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ABSTRACT 

In supernova remnants, the nonlinear amplification of magnetic fields up- 
stream of collisionless shocks is essential for the acceleration of cosmic rays to 
the energy of the "knee" at 10 15,5 eV. A nonresonant instability driven by the 
cosmic ray current is thought to be responsible for this effect. We perform two- 
dimensional, particle-in-cell simulations of this instability. We observe an initial 
growth of circularly polarized non-propagating magnetic waves as predicted in 
linear theory. It is demonstrated that in some cases the magnetic energy density 
in the growing waves, can grow to at least 10 times its initial value. We find no 
evidence of competing modes, nor of significant modification by thermal effects. 
At late times we observe saturation of the instability in the simulation, but the 
mechanism responsible is an artefact of the periodic boundary conditions and 
has no counterpart in the supernova-shock scenario. 

Subject headings: supernova remnants - shock waves - plasmas - cosmic rays 



1. Introduction 



Diffusive shock acceleration (DSA) at supernova remnant shocks is widely co nsidered to 
be the primary source of galactic cosmic rays (for a recent review see lHillasll2005l ). A crucial 
aspect of DSA is the self-excitation of hydromagnetic waves due to currents produced by the 
streaming energetic ions. These provide the pitch-angle scattering necessary to allow spatial 
diffusion of the relativistic particles. The maximum acceleration rate that can be obtained 
corresponds to the smallest possible spatial diffusion coefficient. In th is 'Bohm limit' the 
maxi mum particle energy is determined ei ther by the age of the remnant (ILagage fc Cesarsky 
19831 ). or by its geometry (lBerezhkdll996l ) . Even adopting optimistic values for the upstream 
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parameters of a supernova remnant shock, particles are limited to energies below the knee of 
the cosmic-ray spectrum at ~ 10 15 ' 5 eV. However, nonlinear amplification of the upstream 
magnetic field due to cosmic-ray curren ts may provide a possible solution to this problem 
( Ishapiro et a"Dll998l : iLucek fc Bellhoooh . 



There is increasing observational evidence that the magnetic fields immediately down- 
stream of several young supernova remnant shocks are much stronger t han would be obtained 



by compressing the ambient interstellar field at an MHD shock f ront (IVink fc Laming! 12003 



Berezhko et al. 



2003: 


Bamba et al. 


2005 





additional motivation to the development of t heoretical models of nonlinear field generation 
by both plasma instabilities in the prec ursor jLucek fc Belli l2000l: liell l2004h and fluid-type 



instabilities in the downstream plasma (iGiacalone fc Jokipiil 12007). Th ese two mechanisms 



can, in principle, operate simultaneously (jZirakashvili fc Ptuskin 



20081 ). However, to reach 



the highest energies, particles must be speedily returned to the shock from excursions into 
both the upstream and downstream plasmas. The generation of an amplified magnetic field 
in the precursor that is subsequently advected into the downstream region is sufficient to 
ensure this, whereas amplification of only the downstream field is not. 



Using a hybrid kinetic-MHD analysis, it was demonstrated by iBelll ( 120041 ) that efficient 
cosmic-ray acceleration can, in the linear regime, drive a short wavelength, almost aperi- 
odic, instability with wave-vector parallel to the ambient field. The growth of this mode 
can b e considerably more rapid than that of the resonantly driven modes previously consid- 
ered ( IMcKenzie &: Volk Ill982l ; lAchterbergl Il983l ) . The instability is driven by the reaction 
of the thermal plasma, as it attempts to compensate the current produced by the stream- 



ing cosmic rays. Numerica 



both MHD flBellll2004i 



(PIC) simulations (INiemiec et al 



2005 



investigations of this instability have b een performed using 



Revi 



le et al. 



2008 



2008: 



Zirakashvili et a" 



Riquelme fc Spitkovsky 



20081) and Particle-In-Cell 
20081). While the results of 



MHD simulations show that the perturbed magnetic field S B becomes much lar ger than 
the initial uniform magnetic field Bq, the PIC simulations of INiemiec et all (120081 ) did not 
observe the parallel mode predicted in the linear analysis, and only moderate amplifica- 
tion of the magnetic fie l d was observed. On the other hand, the recent results reported by 
Riquelme &: Spitkovsky! (120081 ) do identify the parallel mode and propose an explanation of 
its saturation in the non-linear stage. 

In this paper we report new PIC sim ulations of th i s inst ability. We extend the range 
of parameters beyond that investigated by INiemiec et all (120081 ). and select values that more 
closely represent the environment in an SNR precursor. Al though both 2D and 3D simu- 
lations in a similar parameter range have been performed (IRiquelme fc Spitkovsky! 120081 ) . 
we restrict our work to 2D. Our results — obtained independently and with some technical 
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differences in the numerical treatment — can be regarded as complementary, and our phys- 
ical interpretation, at least of the early stages of the non-linear development, are congruent. 
In section 2, we recall the linear kinetic analysis of the nonresonant instability, in order to 
elucidate the conditions that must be satisfied in the PIC simulations. The details of our 
simulation parameters are described in section 3. We present the results of the simulations 
in section 4 followed by a discussion of the saturation and its implications in section 5. We 
conclude with a summary of the results. 



2. Linear analysis 



We review the linear kinetic analysis of an electron-ion plasma upstream of a non- 
relativistic shock with a power-law distribution of streaming cosmic rays. We focus our 
attention on parallel shocks, with the streaming velocity directed along the zeroth order 
magnetic field. The kinetic linear dispersion relation for circularly polarized hydromagnetic 
waves prop a gatin g parallel to the mean magnetic field in such a plasma was investigated by 
Achterbergi (119831 ). subject to the conditions of zero net charge and current. In the upstream 
plasma frame the cosmic ray streaming velocity is approximately that of the shock velocity 
( IMcClements et al.lll996l ). Using the diffusion approximation it can be shown that, for a 
power-law momentum distribution of streaming cosmic rays with a spectr al index s > 3 
(fcr(p) oc p~ s ), the dispersion relation can be written as (jReville et al.ll2007l ): 



uT + e 



2 7 2 

oo — v A k 



e(v 2 - (a(kr s ) - 1) = 0, 



(i; 



where u = uo r + vy is the frequency of the waves in the upstream ion rest frame, k > 
the wavenumber, v& the Alfven velocity, and u ci the ion thermal velocity and cyclotron 
frequency, v s the speed of the shock and r g = p m i n c/eB Q the gyroradius of the lowest energy 
cosmic ray. The polarization of the waves is determined by e = 1 ) for right (left) 

handed modes (u r > 0). This result agrees with that found by iBelll (12004 ) using an MHD 



approach, and extends it by including an additional term representing thermal ion damping. 
The growth rate of the unstable modes is determined by the driving parameter 



c 



riim\V s 



(2) 



where n CTj i are the densities of cosmic rays and ions. Finally a is a complex function describing 
the electric fields produced by the cosm ic ray current . It ca n be shown that a is a decreasing 
function of wavenumber for kr g > 1 (jReville et al.l 120071 ). In particular, for wavelengths 
much less than the gyroradius of the cosmic rays, a can be neglected with respect to unity, 
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and one finds that there exists a purely growing mode with growth rate 

C v s v s lv s n cv 
7nr = ~ £ld- (3) 

2 V a r g 2 V a Tl\ 

An important constraint in deriving Eq. (pQ), is the magnetization condition, <C fi c ;. 
On length scales much shorter than their gyroradius, the cosmic rays are essentially un- 
magnetized, i.e., their trajectories are rectilinear. The fact that the background plasma is 
magnetized leads to an asymmetry in the system, resulting in an uncompensated perpendic- 



(Bell 


2005: 


Melrose 


2005) 



support the growing waves, and the ions quickly act to compensate the cosmic-ray current. 
Significant heating of the plasma causes the ions to behave as if they were unmagnetized, 
thereby re ducing the growth-r ate. The same effect arises in the relativistic analysis of this 



instability flReville et al.ll2006h . 



In order to maintain a steady galactic cosmic-ray energy density, supernova remnants 
must channel approximately 10 percent of the kinetic energy entering the shock front into 
cosmic rays. Assuming a shock velocity of v s = 0.01c with a cosmic-ray distribution f CT (E) oc 
E~ 2 from 10 7 eV to 10 15 eV, this implies a cosmic-ray number density of n cr = 10~ 5 cm~ 3 . 
Taking fiducial interstellar medium values for the ion number density and magnetic field: 
m — 1 cm -3 and B — 3 fiG, the maximum growth rate of the nonresonant instability is 
7nr ~ 2.5 x 10~ 3 fi C i. In the upstream rest frame, the cosmic rays drift with respect to the 
thermal ions at the shock velocity v s . To compensate this current, there is relative velocity 
between upstream ions and electrons, Vd = v s n cr /n c = 10~ 7 c. This value is much smaller 
than the ion and electron thermal velocity. Hence an MHD description of the plasma is 
appropriate for the parameters in the precursor of a supernova remnant. P.I.C. simulations 
of this essentially MHD instability are computationally very intensive. However, they are the 
only way of quantifying the relative importance of high frequency waves and other kinetic 
effects on the nonlinear evolution of the instability. 



3. Simulation 



We use a 2D fully relativistic electromagnetic PIC code with a fa st algorithm that 
solves for the current density and conserves charge (jUmeda et al. 1120031). The simu l ations 
are carried out in the upstream rest frame. While the simulations of iNiemiec et all (120081 ) 
used a full kinetic treatment of the cosmic rays, during the linear stage of growth there was 
no significant reduction in the streaming energy of the cosmic rays. In order to identify the 
growth and evolution of the instability, we maintain a constant, uniform external cosmic-ray 
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current throughout the simulations, as has been used in previous MHD simulations (IBell 



2004 ; iReville et all 120081 ). 



Similar to previous simulations, we investigate the role played by streaming cosmic rays 
of uniform density in a homogeneous electron-ion plasma. This is intended to represent a 
small region in the precursor of a supernova remnant shock. However, when, during the 
course of the simulation, the magnetic field fluctuations associated with the nonresonant 
instability become strong, they begin to damp the relative drift velocities of the cosmic rays 
and the gas. In a shock precursor, a change in relative drift velocity is automatically associ- 
ated with an increase in the gas and cosmic ray densities, as the shock is approached. The 
simulation box, however, has periodic boundary conditions and, therefore, simulates a spa- 
tially uniform plasma. These boundary conditions therefore fail to reproduce the behavior 
expected in the shock precursor, as soon as the relative drift velocity changes significantly. 
For this reason, the saturation of the instability observed in the simulation, which occurs 
when the drift speeds become equal, does not correspond to the saturation mechanism ex- 
pected in the shock precursor. However, the initial linear phase of the instability and a 
substantial part of its nonlinear evolution are accurately modelled. 



3.1. Setting 

We define the zeroth order magnetic field to be along the negative x-direction, and 
the cosmic rays stream in the positive x-direction. Periodic boundary conditions are used 
in both the x and y directions. The electrons, ions and cosmic rays are initialized such 
that the overall charge and current densities vanish. Thus, n e = n- x + n CT and the electrons 
have a net drift velocity in the x-direction of V^ e = V^ CT n cx /n e . (Subscripts e, i and cr 
represent electrons, ions and cosmic rays, respectively.) At the beginning of the simulation, 
each population is uniformly distributed in the x — y plane with a Maxwellian momentum 
distribution at equal temperatures T = T e = T; = 1.3 keV. The length of each cell Ax = Ay 
and time step At of the simulation are twice the Debye length and 0.07140;"^, respectively. 
Initially, each cell contains 49 positively and 49 negatively charged macroparticles. In the 
presence of a uniform cosmic ray density, the charge of each species of macroparticle must 
be chosen such that the overall plasma is initially charge neutral. The charge to mass ratio 
of the negatively charged macroparticles corresponds to that of an electron. 

The values for Q cc /u pe , n cr /ni, Vd/c and the number of cells for each simulation run are 
given in Table (U where Q ce and u pe are the electron cyclotron frequency and electron plasma 
frequency, respectively. We also specify the mass ratio m\/m e for each run. To make the 
problem more tractable, we are forced to use unrealistically high values for the cosmic-ray 
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number densities. This is done purely to reduce the computation time, and, as we show, the 
essential physical mechanisms are still well-captured. 

The size of the simulation box in the x and y-directions for run A is taken to be 
L x = L y = 6.66Anr, where Anr = 27ti;a/7nr x i s the wavelength of the most unstable mode 
of the nonresonant instability. The simulation is followed until 40r grow , where r gr0 w = 7nr- 
The parameters used in run B are chosen such that the magnetization condition 7nr < fi C i 
is not satisfied. As discussed in section [2J it is not expected that the nonresonant current- 
driven instab i lity w ill operate in such a situation. Similar parameters have been used by 



Niemiec et all (120081 ). and we perform this simulation in order to facilitate comparisons with 
that work. The size of the simulation box for run-B is L x = L y = 9.65Anr- 

We also perform two additional simulations, runs C and D, with higher cosmic-ray 
drift velocity and mass ratio, respectively. The box size is the same as that used in run A, 
L x = L y = 6.66Anr- While the conditions in run C are less appropriate for a supernova 
shock precursor than those of run A, the increased cosmic-ray drift speed provides a greater 
amount of free energy in the system. Likewise, the magnetization condition in run A is 
satisfied better than in run D, however, the larger mass ratio allows us to investigate the 
effect that increased inertia has on the thermal ions. 



4. Results 

The evolution of the spatially averaged energy density of the magnetic and electric 
field components associated with the growing waves, for run-A is shown in Figure 1. The 
linear growth rate of the spatially averaged transverse component of the magnetic field is 
approximately ~ 0.57nr. The characteristic wave length is almost the same as that predicted 
by the linear analysis. Figure 2a shows the spatial distribution of the z— component of 
magnetic field during the linear phase of growth (t = 14r grow )0 , where the /c-vector is clearly 
directed along the zeroth order field. The y-component of the magnetic field has a similar 
structure with the phase of the wave pattern shifted by tt/2. 

For comparison, we show in Figure 2b, the z— component of the magnetic field during 
the linear phase of growth for run-B. The ions in this case are not magnetized, and the 
nonresonant mode is not observed. The wave is aperiodic and its wavevector is almost 
perpendicular to the beam. This is typical of Weibel-type instabilities. We discuss this 
further in the next section. 



Note t > r grow does not mean nonlinear phase because the initial perturbation is very small. 
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The linear developme nt of the fluid quantities in the simul a tions is similar to that o f 



previous MHD simulations [Belj (12004 120051 ): iReville et all (120081 ): IZirakashvili et all (120081 ) . 



The initial stage sees the development of the aperiodic instability, with uniform density. 
As the fastest growing mode emerges from the initial noise, the net j x B force begins to 
push the plasma in the direction transverse that of the cosmic ray drift, generating low 
density regions between filaments of compressed plasma. The filaments are uncorrelated in 



Bell 


2004; 


Reville et al. 


2008) 



expansion of neighbouring cavities, at which point they appear to have a radius of Anr- 

The late-time nonlinear evolution, however, differs from that of previous MHD simula- 
tions. While both see the eventual disruption of the filamentary structures, the particle in 
cell simulations see an acceleration of the background plasma in the direction of the cosmic- 
ray drift. This effect is not observed in simulations that represent the system as a single 
MHD fluid driven by a fixed external current, because the charge of the streaming particles, 
and hence, the direction in which they drift, is undefined. Figures 3(a) and 3(b) show the 
distribution of the total magnetic field energy and upstream proton density at the end of the 
simulation (t = 40r grow ). Locally, the amplification factor 5B/B is about 16 but the global, 
spatially averaged value is about 4. These images di ffer quite considerably from the no nlin- 
ear results presented in previous MHD simulations (jBell 2004 ; Zirakashvili et al. 2008 ). not 
only because of the net drift, but also due to the anti-correlation between the magnetic field 
and the plasma density (Figure 3(b)). 

At the point at which the magnetic field ceases to grow, the bulk plasma is drifting with a 
constant velocity Vd,i = 0.9vd, cr = 0.09c. This is almost the shock velocity. The development 
of run C and run D are similar to that of run A. However, the saturated field values are 
slightly larger, since the initial amount of free energy available is larger in both cases. The 
growth of the different field components is plotted for run C and run D in Figures 4 and 5, 
respectively. We discuss the physical nature of the field saturation in the next section. 



5. Discussion 

We find that the nonresonant cosmic-ray driven instability develops in all cases where 
the background plasma is magnetized. From Figure [2](a) and [2(b), it can be seen that there 
are considerable differences between the results of run A and run B although the observed 
wavelength is similar to Anr. The most unstable mode in run B is the Weibel type instability, 
occurring between the counterstreaming electrons and ions. Neglecting thermal effects and 
the magnetic field, the growth rate 7^i x an d the wavelength Awi of the Weibel instability 
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where the direction of the wave vector is perpendicular to the drift direction. Comparing 
the two theoretical growth rates and the wavelengths we see that 



.max 

7nr 



— ~ 1, 
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For the choice of parameters used in run B, the wavelengths of the different instabilities are 
qu ite similar, even thou gh the mechanism is quite different. This is consistent with the results 
of iNiemiec et all (120081 ). and we also find the fastest growing mode to be slightly oblique. 
Hence we may be observing a mixed-mode between the Weibel and the nonresonant current 
driven instability or some other type of mixed-mode instability similar to those discussed in 
Bret et al.l J2005h . 



In the late stages of the simulations we see the emergence of a non-thermal power law 
tail in the ion distribution. The use of a constant external cosmic-ray current prevents the 
development of the Buneman instability, a source of ion heating in the linear development. 
T his dramatical l y redu ces the heating in the initial stages, that was previously observed 
in INiemiec et all (120081 ). Figure 6 shows the final electron and ion energy distribution for 



run A, C and D. The distributions are essentially isotropic, and since the bulk drift velocity 
is small in comparision with the thermal velocities, we can safely calculate the spectrum 
in the box frame. The electron energy distribution can be well fitted with a Maxwellian 
distribution with temperature 10 keV for run A and 40 keV for run C and D. A non-thermal 
tail appears in the distribution of the ions with a very soft power-law index of ~ 7.3 for Run- 
A, but a relatively hard ~ 2 for run C and D. The mechanism responsible for the formation 
of this power-law distribution is not clear, and although our simulations have a relatively 
large number of particles per cell, we cannot rule out a numerical artefact. As regions of 
oppositely polarized magnetic field collide each other, an anti-parallel configuration of the 
magnetic fieldlines appears. It is, therefore, possible that collisionless magnetic dissipation 
may also play an important role in plasma heating mechanism and the saturation of the field 
growth as well as the reduction in the drift velocity between the electrons and ions. This 
is an interesting process in its own right, and further investigation is warranted. However, 
in the current context, it is not clear what role these mechanisms play, since the late-stage 
nonlinear behavior has no counterpart in the supernova-shock scenario. 



In agreement with INiemiec et all (120081 ) and iRiquelme fc Spitkovskyl (120081 ) , we observe 
in our simulations that the magnetic field growth ceases when the plasma bulk velocity 
becomes comparable to the cosmic-ray drift velocity. Although we have fixed the cosmic 
ray current in our simulations, the observed saturation is clearly due to the reduction in 



- 9- 



the relative speeds of plasma and cosmic rays. Indeed, if the cosmic rays are isotropic in 
the plasma frame, no streaming instability can operate. The simulations are performed in 
a box with periodic boundary conditions in space. The plasma speed is initially zero, and 
the cosmic ray streaming speed is maintained constant in time and space. It represents the 
speed with which the shock front approaches the simulation box. Within this picture, the 
instability can saturate only when the entire plasma moves with the drift speed of the cosmic 
rays. 

In the precursor of a shock front, conditions are different. There, the interaction between 
plasma and cosmic rays defines the diffusion length scale of the cosmic rays. This is the scale 
on which the cosmic ray current, or, equivalently, pressure, decays with distance ahead of 
the shock front. In the case of efficient acceleration, it is also the length scale on which the 
precursor plasma is compressed. Simulations can model this situation provided the box size 
is small compared to the diffusion length scale, and the time scales are short compared to the 
time on which the cosmic ray current and plasma density change because of the approaching 
shock front. However, this restriction means that the simulations no longer model a precursor 
when they are followed until the plasma speed approaches the shock speed (which equals the 
cosmic-ray drift speed). 

This can be seen explicitly from the equation of conservation of mass in a stationary 
precursor, p(x)u(x) = constant, which dictates that any change in the flow velocity must be 
associated with a change in the density of the fluid. In the present context, the simulation 
box is in the upstream plasma frame, approaching the shock with velocity v s . Defining the 
density and bulk velocity inside the box to be and u\, respectively, we see that 

Pb = 6 

v s - u h 

where po is the initial density of the thermal plasma. Thus, as the flow speed inside the 
box increases, the density should also increase. Since the total density in the box is fixed 
throughout the simulation, the observed growth rate is only accurate provided the bulk 
velocity inside the box is small compared to the shock velocity, or, equivalently, the streaming 
speed of the cosmic rays, i.e., u h <C v s . This corresponds to r grow < 20 for run A and 
r grow < 23 for runs C and D. The amplification of the energy density in the total field in 
run A is not yet substantial, but in runs B and C roughly 10 times the initial energy density is 
reached. The final saturation levels are much higher, but our simulations do not necessarily 
imply that these can be reached in a precursor scenario, and a PIC simulation of this case 
is currently out of reach. 
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6. Summary 

We have performed 2D PIC simulation to investigate the nonlinear physics of the non- 
resonant current instability. We have demonstrated that the energy in the growing waves can 
substantially exceed that of the initial seed field. However, limitations inherent in the simu- 
lation method lead to an artificial saturation level. Although we cannot compute a realisitic 
saturated value of the magnetic field, substantial amplification is seen in some runs before 
they lose validity, and the possibility remains open that the field would continue to grow, 
transferring the energy to longer lengthscales, in accordance with theoretical predictions. 
Thus, as predicted, magnetic field amplification is triggered by the nonresonant streaming 
instability and is likely to play an important role in the acceleration of cosmic rays to the 
knee and beyond in supernova remnants. 
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CfCA, of the National Astronomical Observatory of Japan and on the Blue Gene/P at the 
Rechenzentrum Garching of the Max Planck Society. 
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Fig. I. — Time development of the spatially averaged root mean square of (a) magnetic 
field, (b) electric field for Run-A. Solid, dashed and dotted curves represent the x, y and 
^-components, respectively. 
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Fig. 2. — z-component of magnetic field at (a) t = 14r grow in run-A, (b) t = 8r grow in run-B. 
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Fig. 3. — (a) The magnetic field energy density, and (b) density of ions at t — 40r grow in 
run-A 
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Fig. 4. — Time development of the spatially averaged root mean square of (a) magnetic 
field, and (b) electric field for Run-C. Solid, dashed and dotted curves represent the x, y and 
^-components, respectively. 
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Fig. 5. — Time development of the spatially averaged root mean square of (a) magnetic 
field, and (b) electric field for Run-D. Solid, dashed and dotted curves represent the x, y and 
^-components, respectively. 
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Fig. 6. — Energy distribution at the final stage in run-A, C and D 
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Table 1. Simulation Parameters 



Run 


^ce / ^pe 




m-Jm e 


n CT /ni 


V d ,cr/C 


cell 


A 


3.26 x 10" 


-2 


10 


1/20 


0.1 


4096 x 4096 


B 


4.52 x 10" 


-2 


10 


1/3 


0.3 


512 x 512 


C 


3.26 x 10' 


-2 


10 


1/20 


0.2 


2048 x 2048 


D 


3.26 x 10- 


-2 


40 


1/20 


0.1 


4096 x 4096 



Note. - Parameters of the simulation runs described in 
this paper. Listed are: the ratio of electron cyclotron fre- 
quency and electron plasma frequency Q cc /u pe , ion-electron 
mass ratio m\/m e) the density ratio of ambient ions and cos- 
mic rays n cr /ni, the velocity ratio of cosmic-ray drift and light 
speed Vd, CT /c, the cell number 



